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In this paper we introduce a novel method to simulate lateral diffusion of inclusions in a fluctuating membrane. 
The regarded systems are governed by two dynamic processes: the height fluctuations of the membrane and the 
diffusion of the inclusion along the membrane. While membrane fluctuations can be expressed in terms of 
a dynamic equation which follows from the Helfrich Hamiltonian, the dynamics of the diffusing particle is 
described by a Langevin or Smoluchowski equation. In the latter equations, the curvature of the surface needs 
to be accounted for, which makes particle diffusion a function of membrane fluctuations. In our scheme these 
coupled dynamic equations, the membrane equation and the Langevin equation for the particle, are numerically 
integrated to simulate diffusion in a membrane. The simulations are used to study the ratio of the diffusion 
coefficient projected on a flat plane and the intramembrane diffusion coefficient for the case of free diffusion. 
We compare our results with recent analytical results that employ a preaveraging approximation and analyze 
the validity of this approximation. A detailed simulation study of the relevant correlation functions reveals a 
surprisingly large range where the approximation is applicable. 
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I. INTRODUCTION 

During the last decade it has become apparent that lateral 
diffusion of various components of the cell along the mem- 
brane is crucial for several cellular processes, like exo- and en- 
docytosis, cell signaling, or cell movement. To study all these 
different aspects of lateral diffusion in more detail a whole 
variety of experimental methods has been developed, includ- 
ing fluorescence recovery after photobleaching (FRAP) [1, 2], 
single particle tracking (SPT) [3], fluorescence correlation 
spectroscopy (FCS) |4|], or pulsed field gradient nuclear mag- 
netic resonance (pfg-NMR) |5|]. While the accuracy of mea- 
sured diffusion coefficients achieved with these experimental 
methods can be very high [6], the interpretation of the results 
is often very difficult. It is, therefore, likely that theoretical 
calculations and simulations in particular will play a key role 
in developing a better understanding of diffusive processes in 
biological membranes. 

In order to correctly interpret experimental results it is nec- 
essary to analyze what information of the system in which 
diffusion takes place is documented but also what is neglected 
or insufficiently regarded during the measurement. Lateral 
diffusion of proteins in a membrane must not be viewed as 
diffusion on a flat surface: due to the flexibility of biologi- 
cal membranes lateral diffusion always takes place on curved 
surfaces, whereby the shape of this surface may also be time 
dependent. This is an important aspect that needs to be taken 
into account in the experimental data analysis. However, this 
turns out to be rather difficult, because not always information 
on the membrane shape can be acquired. So far several the- 
oretical studies of free diffusion on temporally fixed curved 
surfaces have been undertaken. One of the first studies of 
this kind was performed by Aizenbud and Gershon J7[] who 
numerically solved the Smoluchowski equation of free diffu- 
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sion on a periodic surface. In later studies that include both 
sophisticated numerical and analytical calculations diffusion 
on more complicated surfaces is reg arded @, & MM O - 
Recently Schwartz et al. lfl3ll and Sbalzarini et al. Qji simu- 
lated free diffusion on surfaces reconstructed from experimen- 
tal image data. In the latter work, that analyses diffusion in the 
endoplasmic reticulum (ER), experimental FRAP curves are 
compared with simulation results. Good agreement is found, 
however, the authors point out that the evaluation of the FRAP 
curves assuming the membrane to be flat leads to diffusion co- 
efficients that differ by a factor of about two from the actual 
intramembrane diffusion coefficient. Furthermore, the experi- 
mentally determined diffusion coefficient appears anisotropic 
while real diffusion along the ER is purely isotropic. 

These studies clearly show that the curvature of membranes 
may not be neglected during the evaluation of experimental 
data. But even if the shape of the membrane within the mea- 
surement volume appears flat on average, the actual shape 
fluctuates due to thermal activation. These membrane fluc- 
tuations have been thoroughly studied. While the fluctuation 
spectrum of a free and almost flat membrane is easily cal- 
culated lfL5i [Trill , the influence of various geometric confine- 
ments of the membrane have also been regarded. These in- 
clude the attachment of the membrane to a reference plane via 
a regular mesh of harmonic springs Odd] which resem- 
bles a model for the attachment of the cell membrane to the cy- 
toskeleton, a membrane that is close to a flat non-impenetrable 
surface 112011 . or the inclusion of active proteins in the mem- 
brane, that exert an out-of -plane force on the membrane l2lll . 
Lin and Brown introduced a very powerful method to simulate 
membrane fluctuations lfl7l [l8l 12211 . the Fourier space Brow- 
nian dynamics algorithm, that allows to add a whole variety 
of external influences on the membrane. Part of the simula- 
tion algorithm introduced in this paper is closely related to 
this method. 

Returning to the movement of proteins along the mem- 
brane, it is obvious that membrane fluctuations will influence 
the measured value of diffusion coefficients. This was first 
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pointed out by Gustafsson and Halle B23L 12411 . But, although 
their work is almost ten years old, experimentally diffusion 
coefficients are still determined by use of the projected flat 
path particles perform instead of the real path along the mem- 
brane. Only recently the quantitative influence of neglecting 
thermal membrane fluctuations on the measured diffusion co- 
efficients has been estimated both by us and Gov ll25l l26ll . 
However, these calculations make use of a preaveraging ap- 
proximation, that a priori implies that the time it takes a diffus- 
ing particle on average to diffuse the length £, should be much 
smaller than the correlation time of a membrane fluctuation 
mode with wave number 2w/E,. In this paper we introduce a 
new algorithm that does not depend on this kind of approxi- 
mation, because the diffusive motion of a protein along a fluc- 
tuating membrane is simulated explicitly. Clearly, the system 
dynamics is described by two coupled processes: the fluctua- 
tions of the membrane and the diffusion of the particle. The 
simulation of membrane fluctuations is effectively regarded 
by numerically integrating the equation of motion for a mem- 
brane given in the Monge gauge, discretely in time. At each 
discrete timestep we also update the position of the diffus- 
ing particle. To this end the discrete version of the Langevin 
equation valid for the particle movement is used. Special care 
needs to be taken regarding this Langevin equation because 
the movement of the particle depends on the actual shape of 
the membrane. 

As a first test we analyze the height-height correlation func- 
tion of membrane fluctuations, and compare the simulation re- 
sults with the analytical result. As mentioned above we have 
previously calculated the ratio of the measured, or projected, 
and the intramembrane diffusion coefficient as a function of 
membrane parameters within a preaveraging approximation. 

Using our simulation scheme that does not rely on this kind 
of approximation we determine the same ratio of diffusion 
coefficients by analyzing the mean square displacement of 
diffusing particles. We compare the simulation results with 
the analytical results and finally discuss the applicability of 
the preaveraging approximation for situations when diffusive 
timescales are comparable or smaller than typical membrane 
time scales by analyzing relevant correlation functions. 

The paper is organized as follows: While we introduce the 
dynamics of membranes in the next section, we explain diffu- 
sion on a curved surface in sec. [Till Generally diffusion can 
either be expressed by a Fokker-Planck equation that gives the 
dynamics of the probability of finding a particle at a certain 
time and position, or a Langevin equation that corresponds to 
the equation of motion of the particle position. In secs. HII A1 
and IIII B I the Smoluchowski equation, a particular form of the 
Fokker-Planck equation, and the Langevin equation for dif- 
fusion on a curved surface expressed in the Monge gauge are 
described, respectively. Since simulation results are compared 
with previous calculations we briefly explain these in sec.lIVI 
After establishing the theoretical foundations necessary for 
this study we describe our simulation scheme in sec. [VI with 
which the results in sec. [Vl] are achieved. In the result sec- 
tion we first analyze pure membrane fluctuations and then de- 
termine the ratio of projected and intramembrane diffusion 
as a function of the membrane parameters bending rigidity 



and effective surface and as a function of the intramembrane 
diffusion coefficient. After the comparison of simulation re- 
sults with analytical calculations we estimate the limits of the 
preaveraging approximation in sec. I VI CI The paper finishes 
with some conclusions and an outlook for future work. 



II. DYNAMICAL MEMBRANE FLUCTUATIONS 

The shape of a membrane without spontaneous curvature, 
i.e., a membrane that is on average flat and without overhangs, 
is conveniently described in the Monge gauge. Hereby a po- 
sition r* of the membrane is given by r* = (x, y, h(x, y)). 
In the following the vector r = (x, y) denotes the projected 
position in the (x, y)-plane. The height function h(r) is the 
distance between the membrane and the flat (x, y) -plane. The 
energy of a membrane with bending rigidity k and effective 
surface tension a is given by 



H[h(r)] = 



r),R 



(1) 



in the Monge gauge 1161 1271 1281 12911 . The projected area of 
the membrane is given by A. A possible energy contribu- 
tion 7Yi[/i(r), R] is caused by the inserted protein at position 
R = (X, Y) and may also depend on the membrane shape. 
In the following we assume that the system's energy does not 
depend on the position of the particle and will therefore drop 
this additional energy term. The dynamics of the membrane is 
given by the equation of motion for the height function h(r). 
Because h(r) is not a conserved order parameter the following 
dynamics applies lfl6l [T7H22I1 : 

= /dV{A(r-r') [ K V 4 Mr',i)-aV 2 Mr',<)]}+£(r,i), 

J A 

(2) 

where A(r — r') is an Onsager coefficient and £(r, t) is a fluc- 
tuating force that obeys the fluctuation-dissipation theorem: 

(CM)) = (3) 

(£(r,i)£(iV)) = 2k B TAA(r-r')S(t-t'). (4) 
Applying the following Fourier transform 



h(k,t) = / d 2 rh(r,t)e 
J A 



A 

k 

the equation of motion for the membrane becomes: 
dh(k,t) 



(5) 
(6) 



Of 



-A(k) [nk 4 + ak 2 ] ft(k, t) + £(k, t), (7) 
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with 

(£(M)) = o (8) 

(£(M)£(k',0) = 2k B TAA(k)6 k ,-v5(t-t'). (9) 

Both the height function h(r, t) and the random force £(r, t) 
are real quantities. Therefore the relation h* (k, t) — h(—\i, t) 
applies (for £ respectively); the asterisk resembles the com- 
plex conjugate. 

The Onsager coefficient for an almost planar membrane can 
be derived from the Oseen tensor and takes the following form 
in Fourier space lfl5l, [l6l, [30ll : 



^ ^ 4:T]k ' 



(10) 



with viscosity r\ of the fluid surrounding the membrane. 

Due to the linearity of eq. (Q the height-height correlation 
function (h(k, t)h(k' , t')), that describes the shape fluctua- 
tions of a membrane is easily calculated fl6ll 

(h(Kt)h(k',t')) = 
nk 3 



exp 



erfc 



4// 



\t-t'\ 



k B TA 
nk 4 + ak 2 



(ID 



III. DIFFUSION ON CURVED SURFACES 

When describing the dynamics of a particle free to diffuse 
laterally along the membrane one has to bear in mind that the 
membrane is not flat but curved. The dynamics of the particle 
is adequately expressed either by a Smoluchowski equation 
that describes the time evolution of the probability of finding 
the particle at a certain position, or by a Langevin equation 
that represents the equation of motion of the particle position. 
In the following two sections both dynamic equations appro- 
priate for a protein diffusing freely on a curved membrane the 
shape of which is given in the Monge gauge, will be intro- 
duced. 



A. Smoluchowski equation 

In Cartesian coordinates the Smoluchowski equation for a 
particle diffusing freely and isotropically on a flat (x, y)-plane 
is given by 



dP'(x,y,t) 
dt 



= DAP'{x,y,t), 



(12) 



with the diffusion coefficient D and the probability 
P'(x,y,t)dx dy of finding the particle in the area element 
dxdy at position (x, y). The probability is normalized such 
that J A dx dy P'(x,y,t) = 1 applies. On a curved surface 
the Laplace operator A needs to be replaced by the Laplace- 
Beltrami operator, which is a function of the metric of the sur- 
face g and the inverse metric tensor g^ |3j]] . The resulting 
Smoluchowski equation takes the form 



dV(x,y,t) 
dt 



1 



Dj^^pd^g^djVix^y^). (13) 



The summations are to be taken over x and y. In the Monge 
gauge the metric of the surface is given by 
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l + ht 



(14) 



with h x = d x h(x, y), for other subscripts accordingly, and 
the inverse metric tensor by 



fJ 



*3 = 



I ' 1 

.9 



hh 



h x hy 
l + hl 



(15) 



The probability V(x, y, t) in eq. ( fT3l ) is normalized such that 
f A dxdyy/gV(x,y,t) = 1 is valid. Assuming that in ex- 
periments typically the path of a particle projected on the flat 
(x, y)-plane is regarded, it makes more sense to evaluate the 
probability P(x, y, t) = -^=V(x, y, t) for which the normal- 
ization J A dx dy P(x,y,t) = 1 applies. To compare results 
that neglect the curvature of the membrane with those that take 
it into account correctly, the differences between P'(x,y,t) 
and P(x, y, t) need to be analyzed. The Smoluchowski equa- 
tion for P(x, y, t) is 



dP(x,y,t) 
dt 



= DY / d i V99 l3 d 3 ^=P(x,y,t) 



(16) 



The probability of finding the projected position of a particle 
within the area element dx dy around position (x, y) is now 
given by P(x, y, t)dx dy. 



B. Langevin equation 

To simulate the movement of a particle in a membrane it is 
more convenient to use a Langevin equation which describes 
diffusion on a curved surface. In general it is possible that sev- 
eral different Langevin equations produce the same dynamics 
of the probability distribution P{x, y, t). In the following we 
will develop a realization of a Langevin equations in the Ito 
calculus that leads to the dynamics of eq. ( fTo*l ). Within the 
Monge description the Langevin equation we wish to develop 
is ideally the equation of motion of the projected position R 
of the particle. The actual particle position is then given by the 
vector (X, Y, h(X, Y, t)). The general form of the Langevin 
equation in Cartesian coordinates is [32] 



dt 



Db, 



GijTj, 



(17) 



where Dbi is the drift term that resembles an external force 
acting on the particle, and GijTj is a stochastic force with 



(Ti(t)) 





26 i:j 6(t-t'). 



(18) 
(19) 



For the curved system we will develop Langevin equations 
that obey the form of eq. (TTTb but where the information on 
the shape of the surface is incorporated into the drift term bi 
and the strength Gij of the stochastic force. 
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The most general form of a Fokker-Planck equation in two- 
dimensional Cartesian space is given by ll32tl 



dP(x,y,t) 

at 



'',!)' 1 0,0 J); P(x,y,t), (20) 



with the drift vector and the diffusive tensor . Com- 
paring this equation with eq. ( [Tol l we can identify 



D™ = D—iWy/gg*)) 
= Dg ij . 



(21) 
(22) 



Note that the partial derivative in eq. (I2T1) is not applied to 
P(x,y,t). If we derive the Langevin equations within the 
Ito calculus the following relationships between the param- 
eters and of the Fokker-Planck equation d20l ) and 
the drift term Dbi and the stre ngth Gu of the stochastic force 
of eq. <H3 need to be fulfilled IBll[32| : 



G, 



D (2 



31 



Dh = D 



(i) 



(23) 
(24) 



Using these relations and the identifications from eqs. (|2TT > and 
d22l we arrive at the following Langevin equation 

ax(t) 

dt 

D-p [2hlh y h xy - h x h xx (l + h 2 ) - h x h yy (l + h 2 x )\ 



+VD——h x hy ( — 



i r„ 



dY{t) 
dt 



D— [2h x hyh xy - h y h xx (l + h 2 y ) - h y h yy (l + h 2 .)] 




i r. 



(26) 



Surprisingly, these equations comprise a drift term that is in- 
duced by the curvature of the membrane and does not appear 
in the Langevin equations for free diffusion on a flat plane. 



IV. FREE DIFFUSION WITHIN THE PREAVERAGING 
APPROXIMATION 

Before we turn to the simulation scheme we will give a 
short introduction to our previous analytical calculations ll25tl 
with which we determine the measured, or projected, diffu- 
sion coefficient of a protein diffusing freely in a membrane. 



The solution of the Smoluchowski equation (fT6b that de- 
scribes the time evolution of the probability distribution of the 
projected position of the diffusing particle, is non trivial be- 
cause the prefactors containing the metric and the inverse met- 
ric tensor are time dependent. If the time it takes a particle 
to diffuse the length £, is much longer than the characteristic 
time r mem b,£, of membrane fluctuations with wavelength £, we 
may apply a preaveraging approximation, i.e., we may replace 
the time dependent prefactors in eq. ( TToT l that contain partial 
derivatives of h(r, t) by their thermal averages. The mem- 
brane time scale is given by the correlation time in eq. ( TTTb 
as r mem b,£ = V^ 3 1 '(27T 3 k), while the diffusive timescale is 
simply given as = E, 2 /4D. Comparing these timescales 
reveals that as long as lengths with 



£,< 7t 3 k/(2D77) 



(27) 



are regarded the preaveraging approximation should be valid. 

When we average over the prefactors of eq. ( fT6] l we find 
that most terms vanish and the Smoluchowski equation sim- 
plifies considerably. Taking into account only leading order 
prefactors it reads: 



dP(x,y,t) 
dt 



= D 



1 



d 2 P(x,y,t) 
dx 2 

hl\ d 2 P(x,y,t) 



a 



dy 2 



(28) 



In an isotropic membrane the two remaining thermal averages 
are both equal and therefore the Smoluchowski equation takes 
the form applicable for diffusion on a flat surface, cf. (fT2l) . 



but now with a new diffusion coefficient D 



proj ■ 



The ratio of 



-Dp ro j that would be measured in experiments, and the actual 
intramembrane coefficient D is given by: 



D 



proj 



D 2 I 1 \ u 



(29) 



By use of the relation 1/g = J* °° da exp[— ag] and the Hel- 
frich Hamiltonian (Q]) this ratio becomes: 



D, 



proj 



D 

1 1 

2 + 2 



da exp 



E 

|k|<q„ 



In 1 



2a 



1 



(3L 2 nk 2 + a 



(30) 



A cutoff wavenumber q m ~ a needs to be introduced that is 
proportional to the inverse of the microscopic length scale a 
in the system. As we will see later, in the simulations this 
microscopic length scale corresponds to the lattice spacing. 
For a more detailed discussion of the numerical evaluation and 
the results of this equation we refer the reader to our previous 
workiH. 

Instead of regarding the preaveraging approximation in 
the Smoluchowski equation it may also be applied in the 
Langevin equations ( f25l ), (l26b . The averaging process leads to 
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a vanishing drift term; the calculation of the the mean square 
displacement and the subsequent derivation of the effective 
diffusion coefficient leads to D ploi /D = {G 2 XX + G^ w ) /2D = 
(1 + (l/g))/2. This is the same result as eq. 



V. SIMULATION SCHEME 

The analytical calculations in the last section were per- 
formed within a preaveraging approximation. In the present 
work we do not apply this approximation, but rather simulate 
the equation of motion for the membrane (O and the Langevin 
equations ( l25T l and ( l26| > for the protein movement. The sim- 
ulation scheme that we will introduce during the next para- 
graphs resembles the effective evaluation of this coupled set 
of dynamic equations by means of discrete numerical integra- 
tion in time. 

First let us turn to the numerical evaluation of the time evo- 
lution of the membrane shape. If we regard eq. (|7]i and as- 
sume periodic boundary conditions it is obvious that the nu- 
merical integration is most effectively implemented in Fourier 
space because the evolution of the height function modes 
h(k, t) takes place independently of each other. If the dif- 
fusion of the protein were not taken into account the whole 
time evolution of the membrane shape could be performed in 
k-space. However, the real space representation h(r,t) be- 
comes necessary to develop the position of the protein. Due 
to the periodic boundary conditions k-vectors are of the form 
k = 2tt(1, m)/L with A ~ L 2 . Because in real space h(r, t) 
is expressed on a quadratic lattice of N x N lattice sites the 
restriction —N/2 < l,m ^ N/2 applies. The lattice spacing 
a is given by a = L/N. 

When numerically implementing eq. (0 one must remem- 
ber that h(k, t) may be a complex number with a real part 
h r (k, t) and an imaginary part hi(k, t) such that h(k, t) = 
h r (k, t)+ihi(k, t). Both for the real and the imaginary part an 
equation of motion is necessary. The numerical equations that 
are used in the simulations to develop the membrane shape 
during a timestep At are of the form 

/i r /i(k, t) = h r/i (k, t - At) 

+ At-^— \nk 4 + ak 2 ] /i r/i (k, t - At) 
Wjk L J ' 



+ y/XksTAAt /(4r)k) r. (31) 

The random number r is Gaussian and therefore (r) = and 
(r 2 ) = 1 applies. The factor A in the random term is either 
1 or 2 as we will now explain: Due to h(r,t) and £(r, t) be- 
ing real quantities not all modes have an imaginary part. In 
particular modes with wave vectors k = 2ir(l,m)/ L with 
(l,m) = (0,0), (0,N/2), (N/2,0), {N/2, N/2) are purely 
real, while all others are complex. In order to fulfill the fluc- 
tuation dissipation theorem from eq. (0 the four purely real 
modes only have an equation of motion for the real part with 
A = 2, while all other independent modes have two equations 
of motion, one for the real and one for the imaginary part, 
each with A = 1. The real and the imaginary part of £(k, t) 
are assumed not to be correlated. Note that not all modes are 



independent, because h(r, t) is a real function. Only a set of 
independent modes is updated via A3 It , while the dependent 
modes are set such that h(k, t) = h*(— k, t). 

Regarding the random term in eq. OTI ) it becomes evident 
that it diverges for k = due to the Onsager coefficient A(k). 
The mode h(k = 0, t) is a measure for the distance between 
the center of mass of the membrane and the flat (x, y)-plane. 
Therefore, fluctuations of h(k = Q,t) just describe a move- 
ment of the membrane as a whole. Such movement of the cen- 
ter of mass is of no relevance for the membrane and diffusive 
properties of interest in this work, so we keep h(k = 0, t) = 
fixed at all times. 

So far the simulation scheme is rather similar to the Fourier 
space Brownian dynamics method introduced earlier by Lin 
and Brown lfl7l fla, l22ll . But in our simulations we addition- 
ally take into account the diffusion of a freely diffusing par- 
ticle along the curved surface given by the membrane shape. 
After an update in the membrane shape using eqs. PTT ) we 
will now update the position of the diffusing particle by using 
a discrete version of eqs. (f25) and (f26) |[3lll32ll : 

X{t + At) = X(t) 

+ [2hlh y h X y ~ h X h XX (l + hl) ~ h X hyy (l + ft*)] At 

D—^—h x hy ( — - 1 1 V2Atr v . (32) 

For Y(t) the corresponding equation is valid. Hereby we use 
h x (R(t), t), for all other partial derivatives of h accordingly, 
at the position of the particle R(t) at time t. The random 
numbers r x and r y are again Gaussian with (r^) =0 and 
( r i r j) = This equation clearly describes an off-lattice 
movement of the particle. It is also conceivable to simulate 
the protein movement through a random walk on the lattice 
used for the membrane. In such a scheme we would check 
the probability of a particle moving to a neighboring lattice 
site within a timestep At and then decide whether the particle 
is to hop. In biological systems, however, the time it takes a 
particle to diffuse the length of a lattice spacing a in the sim- 
ulation is much larger than the typical timescale of membrane 
fluctuations with wave length a. This would typically lead to 
significant changes in the membrane shape before a particle 
jump is successful. Therefore, the off-lattice version of the 
particle's random walk is much more favorable, although the 
partial derivatives of h need to be extrapolated to the position 
of the particle. This is realized by the following procedure: 
Multiplying h(k, t) with the appropriate k-vectors and sub- 
sequently performing the Fourier backtransform, we arrive at 
the first and second partial derivatives of h(r, t) with respect 
to x and y. All discrete Fourier transformations necessary for 
the simulations are effectively implemented using the FFTW 
routines A3 311 . Assume that the protein is to be found some- 
where between the four lattice sites (i, j), (i + l,j), (i,j + l), 
and (i + 1, j + 1), with < i,j < N - 1. The quantity A(R) 
to be extrapolated to the particle position is given at the lattice 
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sites by A(i,j), A(i + l,j), etc. The distance between the 
particle and the line connecting and j + 1) is to be 
/ia and the distance between particle and the line connecting 
(i, j) and (i + 1, j) is va. The linearly extrapolated value of A 
at the particle position R = (a(i + a(j + v)) is calculated 
by use of: 

A(R)=A(i,j) 

- ua(A(i,j) - A(i,j + 1)) - na(A(i,j) - A(i + 
+ ^a 2 (A{i,j) -A(i + l,j) 

+ A(i + lJ + l)-A(i,j + l)) (33) 

After all the necessary quantities have been determined at po- 
sition R the new particle position is calculated with eq. d32i >. 
During the timestep of course also the shape of the membrane 
will change. This is accounted for by employing eq. ( BTT i in 
the next computational step to get h(k, t+ At). Then we again 
update the particle position and so on. Repeating the two 
step process of membrane shape update and particle move- 
ment makes out our simulation scheme. 

Instead of simulating the diffusion of one particle in one 
membrane it possible to insert several particles into the mem- 
brane that do not interact with each other. For each particle a 
separate Langevin equation needs to evaluated, but only one 
membrane equation of motion. Since the Fourier transform of 
the membrane configuration is the most time consuming ele- 
ment of the code, the insertion of more than one particle saves 
computing time. However, the average distance between the 
particles should be sufficiently large in order for the particle 
paths to be independent. Results presented in sections IVIBI 
and lVICl follow from averaging over 2500 paths in 100 differ- 
ent membranes. In each independent membrane 25 particles, 
that do not interact with each other, are allowed to diffuse. 

All simulation results presented in this paper were per- 
formed on a 50 x 50 lattice. If we assume (arbitrarily) that 
the lattice spacing a corresponds to lOnm the system size 
is L = 0.5/im. The viscosity r\ of the water surrounding 
the membrane is r\ = 10~ 9 Js/cm 3 . With the above chosen 
lattice spacing a and the temperature T — 300K, the vis- 
cosity used in the Onsager coefficient of the simulations is 
T) = 2.4 x lQ- 7 k B Ts/a 3 . 

The choice of the appropriate timestep At is determined 
by two timescales, namely the smallest timescale of mem- 
brane fluctuations T mem b,min and the time T<]iff, a it takes a par- 
ticle on average to diffuse a lattice spacing a. Only if At 
is significantly smaller than r mem b,mi n then the evolution of 
the membrane shape will be numerically stable. The mem- 
brane timescale is given by Tmemb.mjn = ^/(nk^ + ak m!lx ), 
cf. eq. ( TTTb . The maximum wave number fc max = \/2n/a is 
determined by the minimal microscopic length scale of the 
system, i.e., in simulations the lattice spacing a. Addition- 
ally the average length a particle moves during At needs to be 
much shorter than a, because we regard the membrane shape 
only right at the beginning of the jump. If the jump is too big, 
the actual particle path along the membrane is not taken into 
account with the necessary accuracy. The diffusive timescale 
is determined by Tdiff. a = a 2 /AD. To perform the simula- 
tions At should be considerably smaller than both T mem b >m in 




k [units of 1/a] 

FIG. 1: Height-height correlation function {h(k)h* (k)} as a func- 
tion of wave number k (in units of 1 /a). To illustrate the agreement 
between the analytical and the simulation result for larger k values, 
we plot k 2 / {h(k)h* (k)) as a function of k 2 in the inset. The sim- 
ulation results, symbolized by the diamond symbols, were obtained 
by averaging over 250 independent membrane configurations cre- 
ated during runs with dimensionless physical parameters f3n — 5 
and (3a L 2 = 500 on a 50x50 lattice. The numerical timestep was 
At = 1.7 x 10~ 10 s. 



and Tdiff.a. Timesteps At used in the presented simulations 
range from 5 x 10 _10 s to 2 x 10 _9 s. Typical simulation runs 
comprise ~ 2 x 10 6 timesteps which took approximately 30- 
35 minutes on a 64 node Beowulf cluster with Pentium IV, 3.2 
GHz processors. 



VI. RESULTS 
A. Validation of membrane fluctuations 

After developing a simulation code it is always necessary 
to test it by comparing its results with results previously ob- 
tained with another method. In this section we will show that 
the evolution of the membrane shape results in the membrane 
fluctuations given by eq. ( fTTT i that follow analytically for a 
membrane with the Helfrich energy (Q]) and the Onsager coef- 
ficient of eq. ( [Tol l. 

In fig. Q] we display the equal time correlation function 
(/i 2 (k, t) + /if(k, f)) as a function of k = |k| for (3k = 5 
and (3aL 2 — 500. This corresponds to a surface tension on 
the order of 8 x 10~ 3 mJ/m 2 . Typical values for the dimen- 
sionless bending rigidity j3n of lipid bilayer membranes are 
between 1 and 50. 

If we compare the simulation results with the analytical re- 
sult that is given by the solid line we see a good agreement 
for small wave numbers. For higher k values it is more con- 
venient to plot k 2 j (hl(k, t) + ft|(k, t)) as a function of k 2 , 
which is done in the inset. The linear behavior expected from 
the analytical calculations is well reproduced by the simula- 
tions. 
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FIG. 2: Time correlation function (h(k, t)h*(k, 0)) as a function of 
the dimensionless time t/r mcm b(fc) for the given values of k. The 
correlation time r mcm b(fc) = 4ry/(Kfc 3 + ak) is given by eq. i ll It , 
Symbols are results from the simulations, while straight lines resem- 
ble fitS OC exp[— t/Tmemb]- 



FIG. 3: (Color online) Elucidation how D pr0 j is determined from sim- 
ulation results (symbols): the averaged mean square displacement 
(AR 2 (t)) /4D is plotted as a function of time. The error bars cor- 
respond to the standard error of averaging. The slopes of linear fits 
(solid lines) to these results determine the diffusion coefficient. The 
projected diffusion coefficient on a fluctuating membrane is always 
smaller than on a flat plane (thick dashed line). 



In order to test the time correlations of the membrane fluc- 
tuations we plot the correlation function (h(k, t)h* (k, 0)) as a 
function of the dimensionless time t/r mem b(fc) for several ran- 
domly chosen wave numbers in fig. [2] The same parameters 
as in fig.[T]were used. The exponential fits oc exp[— i/r mem b] 
to the simulation results reveal a good agreement between the 
analytical correlation times and the results obtained from our 
numerical membrane update scheme. 



B. Free membrane-bound diffusion 

In the previous section we showed that our simulation 
scheme reproduces membrane fluctuations correctly. In this 
section we show that it is also capable of adequately describ- 
ing the diffusion of a protein in the membrane. 

To determine the projected diffusion coefficient in the sim- 
ulations we use the relation ((R(t) - R(0)) 2 ) = 4D prq ji. The 
mean square displacement (AR 2 (t)) = ((R(i)— R(0)) 2 ) as a 
function of time is determined from the simulations by averag- 
ing over a large number of independent particle paths. In fig. [3] 
we show (AR 2 (i))/4Z? as a function of time for the given 
values of (3k and a = 0. The chosen intramembrane diffusion 
coefficient is D = 10 5 a 2 /s; this corresponds to an experi- 
mental value of D = 10~ 7 cm 2 /s. The results evidently show 
a linear increase of (AR 2 (i)) with time. Furthermore, the 
comparison of the simulation results with the line that gives 
the expected mean square displacement if the membrane were 
flat displays that the projected diffusion coefficient is smaller 
than the actual intramembrane diffusion coefficient. This is of 
course expected and also seen in the result of eq. j29) , because 
membrane fluctuations increase the actual path of the protein. 

The slope of the linear fit to the simulation results some of 
which are displayed in fig.[3]corresponds to the ratio D pro j/D. 
In fig. |4]we plot D pm j/D resulting from the simulations and 
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FIG. 4: Comparison of simulation and analytical results for D pm j/D 
as a function of bending rigidity {3k for the two given effective ten- 
sions. Good agreement is observed. 



the numerical evaluation of the integral in eq. d30i > as a func- 
tion of bending rigidity (3k both for vanishing tension and 
(3<jL 2 = 500. Both simulations and preaveraging calcula- 
tions show that for small bending rigidity and small tension 
the difference between projected and actual diffusion coef- 
ficient becomes more and more pronounced. This result is 
plausible: decreasing tension a and bending rigidity k leads 
to stronger membrane fluctuations. The difference between 
actual and projected particle path becomes bigger and there- 
fore the fluctuation effect is enhanced. If the system size is 
increased, not shown here, this also leads to an increase in 
fluctuation strength. The comparison of simulation and ana- 
lytical results displays a very good overall agreement. 

Results in fig. [4] were achieved for a single diffusion coef- 
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FIG. 5: (Color online) Ratio of the projected to the intramembrane 
diffusion coefficient D pm j/D as a function of (3k for various diffu- 
sion coefficients (given in units of a 2 Is) and a = 0. Simulation 
results for all D cannot be distinguished from the preaveraging re- 
sult. 



ficient D. To study the influence of the intramembrane diffu- 
sion coefficient we also perform simulations for five different 
coefficients D = 10 4 , 10 5 , 10 6 , 5 x 10 6 , 10 7 a 2 /s and a = 0. 
In fig. [5] we display the ratio D pm j/D as a function of (3k for 
all used D as derived from the average mean square displace- 
ment (AR 2 (i)) of the diffusing particles. The results show 
that D pm j/D is seemingly independent of the diffusion coef- 
ficient D: for all values of D a good agreement of the sim- 
ulation results is observed. Comparing the simulation results 
with the analytical result invoking the preaveraging approx- 
imation we also find that the agreement is very good. But 
this is rather surprising: a priori the preaveraging approxi- 
mation should only be applicable if the time it takes a par- 
ticle to diffuse the length £, is much larger than the correla- 
tion time of membrane fluctuations with wavelength £,. This 
is expressed in eq. d27l ) where we have a crossover length 
scale f, co ee 7t 3 k/(2D?7) ~ 6 x 10 7 /3n/D. If diffusion on 
length scales below £, co is analyzed the preaveraging approx- 
imation should lead to good agreement, while one would ex- 
pect to see a change for larger length scales. For the small- 
est regarded diffusion coefficients £, co is on the order of 10 3 a. 
With a system length of L = 50a all lengths in the system 
are below the crossover length scale and therefore the agree- 
ment of the simulation results with the preaveraging result, as 
seen in figs. [4] and [5] is expected. For the largest regarded D 
the crossover length is on the order of 10a. Thus there are 
very many membrane fluctuation modes in the system with 
wave lengths larger than £, co . In this case we would have ex- 
pected the interplay of membrane and diffusive timescales to 
be observable in the effective diffusion coefficient. However, 
fig. shows that this is not the case. To understand why the 
preaveraging approximation leads to good agreement with ex- 
plicit simulations even for large D, it is necessary to study the 
correlations of the drift term in the Langevin equations ( f25T > 
and d26*l ) that is caused by the metric of the system. This is the 
subject of the following section where we discuss the validity 



of the preaveraging approximation. 

C. Validity of the preaveraging approximation 

The last section had the unexpected outcome that the ana- 
lytical calculation within the preaveraging approximation de- 
scribes particle diffusion well even when the diffusion coef- 
ficient is so large that one would assume this approximation 
to break down. In this section we study in more detail the 
validity of the preaveraging approximation. 

To this end it is helpful to regard the Langevin equa- 
tions d25l l and d26l ) governing particle diffusion. The diffusion 
coefficient is defined through the mean square displacement 
of the diffusing particle. Using eq. ( fT7b we can formally write 
the mean square displacement as 

(AR 2 W> = 

D 2 J* dr J* dr' j^^RCr); r)6 j (R(r / ); r')) j 

+ 2j2(Gl)t (34) 



The last term that is linear in time t coincides with the preav- 
eraging result, because (G 2 XX + G 2 yy )/2D = (1 + (l/g))/2, 
see eq. d29l ). Hence the diffusion coefficient derived through 
(AR 2 (<))/4i gets an additional term caused by the correla- 
tions of the drift term. If we are interested in the ratio of 
the projected and the actual diffusion coefficient it may be 



written as the sum of two terms D pm j/D = D p 



v/D 



D 



proj, drift 



/ D, with the additional contribution 



D, 



proj, drift 



j t D J' dr ^ dr' {(b x [h(R(r);T)]b x [h(R(r');r')}) + 
{b y [h(R(ry,T)]b y [h(K{T'y,r')})}, (35) 

with 

b x [h} = -^ ■2lrli, i h : „ /',//,, {l + h 2 y )-h x h yy {l + h 2 x )\ , 

b y [h] = -2 [2h 2 h x h xy -h y hyy (l + hl)-h x h xx (l+h y )] . 

(36) 

In order to estimate this contribution we calculate the 
functions (& x (R(r); t)6 2; (R(t + Ar);r + At)) and 
(6 y (R(r); t)6 2/ (R(t + Ar);r + At)) from our simulation 
data for the five previously regarded diffusion coefficients. 
Note that the correlations do not only depend on the pure 
time interval |r — t'\ but also on the distance |R(r) — R(t')| 
the particle travels during this time interval. The correla- 
tion functions are displayed in the top and bottom panels 
of fig. [6] as a function of At for (3k = 1 and f3n = 2. 
In this whole section we set the effective surface tension 
(7 = 0. We find that the correlation function is overall smaller 




FIG. 6: (Color online) Correlation functions {fe a; (R(r); r)6 a; (R(r + 
At); t+At)) (solid symbols) and (6 H (R(r); r)6 y (R(r+ At); t + 
At)) (white symbols) as a function of At for (3k = 1 (top panel) 
and (3k — 2 (bottom panel). Different symbols apply for different 
diffusion coefficients (given in units of a 2 /s). The insets display the 
same results as double logarithmic plots. 



for larger (3k as is expected. This means that the larger 
the bending rigidity n or the effective tension a, not shown 
here, the less important is the contribution caused by the drift 
term. The expectation for an isotropic system that the corre- 
lation functions (6 x (R(t); t)6 k (R(t + At);t + At)) and 
(& y (R(r); r)by(R.(T + At);t + At)) coincide, is also ful- 
filled. Furthermore, a faster decrease of the correlation func- 
tion with increasing diffusion coefficient D is observed. The 
decrease of the correlation function with time At is deter- 
mined by two processes: the change in membrane shape and 
the movement of the particle along the membrane. If the parti- 
cle did not diffuse, i.e., D = 0, the displayed correlation func- 
tions would still decrease with increasing the time interval due 
to the evolution of the membrane shape. This contribution to 
the decrease obviously does not depend on how fast the parti- 
cle diffuses within the membrane, but only depends on mem- 
brane parameters like bending rigidity k or surface tension a. 
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1 2 3 4 5 

DAt [units of a 2 ] 



FIG. 7: (Color online) Correlation functions (b x {R,(r); T)b x (R.(T + 
At); t+At)) (solid symbols) and {b y {R{r); t)6 h (R(t + At); t + 
At)) (white symbols) as a function of DAt and At (inset) resulting 
from diffusion on 100 different quenched membranes with 25 parti- 
cles each for (3 k = 1. Different symbols apply for different diffusion 
coefficients. 



In the other extreme when membrane time scales T mem b ap- 
proach infinity, and the particle diffuses on a fixed membrane 
shape, the decrease of the correlation function is determined 
by the distance a particle travels during a fixed time interval. 
This is shown in fig. [7] for (3k = 1, where we display the 
correlation functions (&i(R(r); t)&j(R(t + At);t + At)) 
for particles diffusing on different quenched membranes over 
which we average afterwards. The "quenched" configurations 
used in the simulations are obtained by evolving the mem- 
brane shape for such a long time that thermal equilibrium has 
been reached. Regarding the correlation functions as a func- 
tion of time we see that the smaller the diffusion coefficient 
the slower the decrease of the correlations. Multiplying At 
with the diffusion coefficient leads to a perfect match of all 
four lines. This is a clear indication that the correla tion func- 
tion depends only on the average distance V AD At a particle 
travels during a certain time. The t hick sol id line is the re- 
sult of an exponential fit oc exp(— V4DAt/C); the correla- 
tion length for (3k = 1 is given by £ ~ 1.0a. 

If we return to fig. |6]where both the membrane and the par- 
ticle are moving, it is now understandable that an increasing 
diffusion coefficient leads to a faster decrease of the correla- 
tion functions. An analytic form for the observed correlation 
functions could not be found. 

In order to estimate the additional contribution to the pro- 
jected diffusion coefficient the correlation functions need to 
be integrated according to eq. ( [34-b . Regarding the results in 
fig- E] it is obvious that the double integration of the correla- 
tion functions (&j(R(t); T)6j(R(T + At); t + At)) in time 
will cause a slower increase with time for large D. However, 
the additional contribution -D pro j, drift /D is given by multiply- 
ing this integral by D. The additional mean square displace- 
ment (AR^ rift (i))/4Z) resulting from the numerical integra- 
tion of the results in fig. [6] is displayed in fig.[8]for (3k = 1. 
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FIG. 8: (Color online) Additional mean square displacement 
(ARjm(t))/4.D as a function of time f for (3k, = 1 as derived by 
numerically evaluating the first term of eq. ( 134-t . 



For all regarded diffusion coefficients we find a linear behav- 
ior for large times t. Furthermore, the slope becomes larger 
for increasing D. In other words the influence of the drift 
term becomes more pronounced the faster the particle diffuses 
compared to the fluctuations of the membrane. 

In this context it is interesting to check whether the addi- 
tional contribution is finite for an infinite intramembrane dif- 
fusion coefficient. For very large diffusion coefficients the 
membrane will appear almost stiff for the moving particle. 
This corresponds to the situation regarded in fig. [TJ where we 
found that the correlation function is well described by an ex- 
ponential function (6 j; (R(r); r)6j(R(r + Ar);r + At)) ~ 
(bf (R(r); r)) exp[-V4X>A~77C]. For this particular function 
eq. ( f35l > is easily evaluated. In the limit of large times t we 
find D ?IOhdm /D = (6f(R(r);r))C 2 /2. Thus D pmiM& /D is 
independent of D, and only involves membrane parameters. 
It therefore remains finite. For (3k, — \ the largest possible 
increase of D pm j/D is approximately ~ 0.015. 

The additional terms to the ratio of projected to intramem- 
brane diffusion coefficient -D pr qj .drift /-D resulting from fits to 
(ARj ritt (i)), see fig. [8] are displayed in fig.|9] Although an in- 
crease in D causes a larger additional term we still see that for 
the regarded diffusion coefficients that are much larger than 
those experimentally observed in experiments, the additional 
term is always more than two orders of magnitudes smaller 
than -D P roj,preav/-D- For D — 10 7 a 2 /s the numerically deter- 
mined value of -Dp r oj,drift/-D agrees reasonably well with the 
previous estimate for infinite D. Surprisingly, for situations 
when one expects the preaveraging approximation to break 
down, it still gives reliable results. Even for infinite diffusion 
coefficients the results for the projected diffusion coefficient 
calculated within the preaveraging approximation only differ 
from the actual values by less than two percent for experimen- 
tally accessible membranes. 

After this discussion we can understand why our simula- 
tion runs agree well with the preaveraging result for all D, as 
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FIG. 9: Additional contribution -D P ioj,drift/-D following from linear 
fits to (ARj rift (f)), see fig. [8] as a function of D for (3k = 1, 2. 



can be seen in fig. [5] While our simulation scheme is able to 
achieve an accuracy of a few percent within reasonable com- 
puting time, the corrections caused by the drift term cannot 
be identified directly via the mean square displacement of the 
particles. Nevertheless, the explicit evolution of the mem- 
brane shape and particle position make it possible to estimate 
this correction via the evaluation of correlation functions. 



VII. CONCLUSIONS 

In this paper we introduce a novel scheme that allows for 
the simultaneous simulation of membrane fluctuations and in- 
tramembrane diffusion. The dynamics of the system is ex- 
pressed via the equation of motion for a membrane described 
in the Monge gauge and the Langevin equation for a particle 
diffusing along a surface whose form is given in the Monge 
gauge. The simulation algorithm consists of the numerical 
integration of these two coupled differential equations. To 
validate the membrane fluctuations we compare the height- 
height correlation function determined from the simulations 
with known analytical results. After ensuring that membrane 
fluctuation are reproduced correctly we study free membrane 
bound diffusion along the membrane. Since diffusion coeffi- 
cients are experimentally often determined from the projected 
path a particle covers, we regard the ratio of the measured, 
projected diffusion coefficient and the intramembrane diffu- 
sion coefficient that is a parameter of the simulations and cal- 
culations. Both the simulations and the previous calculations 
that apply a preaveraging approximation, show that the differ- 
ence between the measured and the true intramembrane dif- 
fusion coefficient is largest for small bending rigidities k and 
small effective surface tensions a. This can be understood be- 
cause small k and a lead to stronger membrane fluctuations. 
Thus the actual path and the projected path differ most. Our 
calculations reveal a maximum reduction of the projected dif- 
fusion coefficient by approximately 20 percent. We are aware 
that the experimental corroboration of our findings is currently 
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challenging, but with the constantly increasing accuracy of 
methods to determine lateral diffusion it should become fea- 
sible in the near future. Such experiments will be important 
in showing that lateral diffusion is not only a function of the 
direct interaction of lipids and proteins but also depends on 
material properties of the membrane. 

We also consider simulation runs with different intramem- 
brane diffusion coefficients D. The subsequent analysis re- 
veals a surprising observation: the resulting ratios D pm j/D all 
coincide independently of D. Furthermore, the simulation re- 
sults agree well with the analytical preaveraging calculations. 
Only simulation runs with the smallest regarded diffusion co- 
efficients D are expected to be well described by the analytical 
results. For the largest used D, however, when diffusive and 
membrane time scales become comparable, one would a pri- 
ori assume that the ratio D pm j/D from the explicit simulations 
would differ from the calculations. 

An analysis of the Langevin equation that determines the 
movement of the particle, demonstrates that correlations of the 
drift term caused by the metric of the membrane are respon- 
sible for a possible increase in the measured diffusion coeffi- 
cient. In order to understand the applicability of the preaver- 
aging approximation we study these correlation functions us- 
ing our simulation scheme. The relevant correlations decrease 
in time not only due to membrane fluctuations but also due 
to the movement of the particle. Our simulations reveal that 
the influence of the drift term on the projected diffusion coef- 
ficient increases with increasing D. But surprisingly for all, 
even infinite, intramembrane diffusion coefficients it is very 
weak in experimentally accessible membranes. In fact the in- 
fluence is so small - in our simulations below two percent -, 
that it cannot be directly identified from studying the mean 
square displacement of the particles within our scheme. Only 
the study of the correlation functions of the drift term using 
our simulations gives insight into the additional contributions 
to the projected diffusion coefficient. For future studies one 



may now argue that preaveraging suffices and the simulation 
scheme becomes unnecessary. Note that this reasoning only 
makes sense as long as the particle diffuses freely. If an addi- 
tional interaction between membrane and protein is included 
the analytical calculation of both the altered diffusion coeffi- 
cient and membrane spectrum relies on approximations which 
the simulations are capable of overcoming. 

The experimental study of lateral diffusion in membranes 
has become a very important and much noticed field that 
has revealed many different diffusion phenomena. But so far 
only little complementing theoretical or simulational work has 
been performed in order to develop a broader understanding 
of experimental findings. In previous work B25I1 we regarded 
the influence of a simple interaction between a diffusing par- 
ticle and the membrane curvature and found a strongly altered 
dependence of the projected diffusion coefficient on bending 
rigidity and effective tension. It is therefore promising that 
the measurement of lateral diffusion coefficients as a function 
of membrane properties might shed light on the interaction be- 
tween protein and membrane. While our previous calculations 
employed several approximations the method introduced in 
this paper overcomes these. Our powerful simulation scheme 
is, therefore, a good starting point to investigate the influence 
of various membrane-protein interactions not only on lateral 
diffusion but also on the spectrum of the membrane. Further 
effects that are easily incorporated into our scheme are ex- 
ternal influences on the membrane, like tethers that resemble 
the attachment of a membrane to the cytoskeleton. These and 
other extensions will be part of our future work. 
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